function IdSet=OneSide(u0grid, u1grid, u2grid, u3grid, u4grid, ...
                       sgreduced, PYX, ...
                       id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
                       P_col,...
                       Tcoord_zeros,...
                       Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
                       cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros, equalorder, n_vertices,infty_pos)

exists=zeros(sgreduced,1);
param_MOSEK.MSK_IPAR_LOG = 0; 
param_MOSEK.MSK_IPAR_INTPNT_BASIS = 'MSK_BI_NEVER';
param_MOSEK.MSK_IPAR_NUM_THREADS = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Outside the loop across parameter values %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Observational equivalence
% RHS 
RHS_obsequiv=(-1+PYX).';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
%%%%%%%%%%%%%%%% Inside the loop across parameter values %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for g=1:sgreduced  

% Set indices to keep track of equal elements     
idU1=id_Ugridreduced{1}(g,:);
idU2=id_Ugridreduced{2}(g,:); 
idU3=id_Ugridreduced{3}(g,:); 
idU4=id_Ugridreduced{4}(g,:); %for each of the U-sets
% Relevant sizes of the vector of unknowns
s1=s_id_gridreduced(g,1);
s2=s_id_gridreduced(g,2);
s3=s_id_gridreduced(g,3);
%s4=s_id_gridreduced(g,4);
sU=sU_id_gridreduced(g); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Equalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Observational equivalence
% Coordinate columns
cc_obsequiv=[changecoord_4(s1,s2,s3, idU1(P_col{1}(:,1)), idU2(P_col{1}(:,2)), idU3(P_col{1}(:,3)), idU4(P_col{1}(:,4))) ... %Y=1
             changecoord_4(s1,s2,s3, idU1(P_col{2}(:,1)), idU2(P_col{2}(:,2)), idU3(P_col{2}(:,3)), idU4(P_col{2}(:,4))) ... %Y=2
             changecoord_4(s1,s2,s3, idU1(P_col{3}(:,1)), idU2(P_col{3}(:,2)), idU3(P_col{3}(:,3)), idU4(P_col{3}(:,4))) ... %Y=3
             changecoord_4(s1,s2,s3, idU1(P_col{4}(:,1)), idU2(P_col{4}(:,2)), idU3(P_col{4}(:,3)), idU4(P_col{4}(:,4))) ... %Y=4
             changecoord_4(s1,s2,s3, idU1(P_col{5}(:,1)), idU2(P_col{5}(:,2)), idU3(P_col{5}(:,3)), idU4(P_col{5}(:,4)))];   %Y=0

%% Integrates to one
cc_one=changecoord_4(s1,s2,s3, idU1(infty_pos), idU2(infty_pos), idU3(infty_pos), idU4(infty_pos));

%% Zeros
cc_zeros=changecoord_4(s1,s2,s3, idU1(Tcoord_zeros(:,1)), idU2(Tcoord_zeros(:,2)), idU3(Tcoord_zeros(:,3)), idU4(Tcoord_zeros(:,4))); %1xn_zeros


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Inequalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Increasingness 
%Repeat the way we constructed Tcoord but now using the parameter values under consideration
vectors = cellfun(@(x) {x(g,:)}, Ugridreduced); %cell 1xn_U_sets
%For j=1,...,n_U_sets, it picks the g-th row of Ugridreduced{j}
T_temp = cell(1,n_U_sets);
[T_temp{:}] = ndgrid(vectors{:}); 
T_temp = cat(n_U_sets+1, T_temp{:});
T = reshape(T_temp,[],n_U_sets); %matrix of dimension as Tcoord reporting the elements of all the possible n_U_sets-tuples from the U-sets

D=[T(indices_pairs_extended(:,1),:) T(indices_pairs_extended(:,2),:) ...
   Tcoord(indices_pairs_extended(:,1),:) Tcoord(indices_pairs_extended(:,2),:)]; %n_U_sets-tuples n_U_sets-tuples coordinates coordinates


%Save in D1 the rows of D where the first n_U_sets-tuple is <= than the second
%Only the columns reporting column-coordinates are relevant
D1=D(sum(D(:,1:n_U_sets)<=D(:,n_U_sets+1:2*n_U_sets),2)==n_U_sets, 2*n_U_sets+1:end); %sd1x(2*n_U_sets)
sd1=size(D1,1);

%Save in D2 the rows of D where the first n_U_sets-tuple is >= than the second
%Only the columns 5:end (reporting column-coordinates) are relevant
D2=D(sum(D(:,1:n_U_sets)>=D(:,n_U_sets+1:2*n_U_sets),2)==n_U_sets, 2*n_U_sets+1:end); %sd2x(2*n_U_sets)
sd2=size(D2,1);

%Flip D2
D2=[D2(:,n_U_sets+1:2*n_U_sets) D2(:,1:n_U_sets)]; %first n_U_sets-tuple is <= than the second

%Combine D1 and D2
D=[D1;D2];
sd=sd1+sd2; 

cr_increas=kron((1:1:sd), ones(1,n_vertices)); %1x[sd*n_vertices]
%n_vertices comes from the number of all possible vertices generated by each pair of n_U_sets-tuples

cc_increas=zeros(1,sd*n_vertices);
for j=1:sd
    D_temp=D(j,:);
    vertices=D_temp(C); %n_verticesxn_U_sets
    cc_increas(n_vertices*(j-1)+1:n_vertices*j)=changecoord_4(s1,s2, s3,idU1(vertices(:,1)), idU2(vertices(:,2)), idU3(vertices(:,3)), idU4(vertices(:,4))); %1xn_vertices
end

fill_increas=repmat(-C_sign, 1, sd);  %1x[sd*n_vertices]
%remember: signs are flipped because inequalities are <= 

RHS_increas=zeros(sd, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Compose %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Aeq=sparse([cr_obsequiv cr_zeros cr_one ], [cc_obsequiv cc_zeros cc_one ], [fill_obsequiv fill_zeros fill_one], numeqconstraints, sU);   %[#equalityconstraints]x[sU]
Aineq=sparse(cr_increas, cc_increas, fill_increas, sd1+sd2, sU);   %[#inequalityconstraints]x[sU]
beq=[RHS_obsequiv; RHS_zeros;RHS_one];

% MOSEK
prob.blx=zeros(sU,1); %lower bound unknowns
prob.ulx=ones(sU,1); %upper bound unknowns
prob.c=zeros(sU,1); %objective function
prob.a=[Aeq; Aineq];
prob.blc=[beq; -Inf(size(Aineq,1),1)];
prob.buc=[beq; RHS_increas];
[~,result]     = mosekopt('minimize echo(0)',prob, param_MOSEK);
if strcmp(result.sol.itr.prosta, 'PRIMAL_AND_DUAL_FEASIBLE') %if the primal is feasible
   exists(g)=1; 
end

end

%Enlarge results to the original grid of parameter vectors
exists=exists(equalorder);
IdSet = [u0grid(logical(exists),:)  u1grid(logical(exists),:)  u2grid(logical(exists),:) u3grid(logical(exists),:) u4grid(logical(exists),:)];

end

